The Optimal Centroidal Voronoi Tessellations and the Gersho's Conjecture in the Three-Dimens ional Space

K e y w o r d s O p t i m a l triangulation, Centroidal Voronoi Tessellation, Gersho's conjecture, Optimal vector quantizer, Mesh generation and optimization. We would like to thank the referee for many helpful suggestions that improved our presentation. The research was supported in part by the China State Major Basic Research Fund G199903280 and by the NSF-DMS 0196522. 0898-1221/05/$ see front matter @ 2005 Elsevier Ltd. All rights reserved. Typeset by Jt3dS-~_tX doi: 10.1016/j .camwa.2004.12.008 1356 Q. Du AND D. WANG 1. I N T R O D U C T I O N A centroidal Voronoi tessellation (CVT) is a Voronoi tessellation of a region such that the generating points of the tessellations are also the mass centroids of the corresponding Voronoi regions with respect to a given density function [1]. Its dual triangulation is called the centroidal Voronoi Delaunay triangulation (CVDT) [2,3]. CVTs enjoy an optimization characterization so that they turn out to be useful in diverse applications such as image and data analysis, vector quantization, resource optimization, optimal placement of sensors and actuators for control, cell biology, territorial behavior of animals, high quality mesh generation, numerical partial differential equations, meshless computing, etc., [1,3-7]. Given a domain, a specific density and a fixed number of generators, the optimal CVT is defined as one which has the lowest cost functional (or energy, distortion) among all the possible CVTs. Such configurations are directly related to the celebrated Cersho conjecture. The conjecture was originally stated for optimal vector quantizers used in data compression and transmission [8], and translating in the language of CVTs [1], it states that asymptotically as the number of generators gets larger and larger, the optimal CVT will be forming a regular tessellation consisting of the replication of a single polytope whose shape depends only on the spatial dimension. The basic Voronoi cell of the optimal CVT was shown to be the regular hexagon in two dimensions [9]. For the three-dimensional space with a constant density function, it was proved that among all lattice-based CVTs, the CVT corresponding to the body-centered cubic lattice (BCC) is the optimal one [10]; while for nonlattice cases or the general cases, the conclusion remains to be rigorously proved [8,11]. In this paper, we conduct some well-organized numerical experiments on CVTs to provide a clear substantiation of the Gersho conjecture in 3D: the BCC lattice based CVT is not only optimal among the lattice based ones, but also is likely the global optimal solution for all the CVTs including nonlattice ones. The computations are done mostly via the classical Lloyd iterations [12]. While the convergence of such iterations to a CVT coniigurration is mostly assured [13], its convergence to the global optimum is dependent on the choice of the initial configuration. We demonstrate that the optimal CVT and its dual CVDT enjoy very nice geometrical and topological properties, and they lead to high quality tetrahedral meshes in the three-dimensional space. For simplicity, much of our simulations are carried out for the constant density case since we are mostly interested in establishing results in the asymptotic regime where the variation of the density in space is less important. Such a simplification also allows us to conduct numerical experiments on relatively smaller scales than those otherwise required for more complicated nonuniform densities. In addition, we also design careful experiments in various geometry to address the boundary effect on the simulation results. The remainder of the paper is organized as follows: first in Section 2, we discuss how to compute the optimal CVTs, mostly based on the Lloyd iterations. In Section 3, abundant examples are presented to show that the BCC based CVT is likely the optimal CVT for both lattice and nonlattice based configurations. In Section 4, we illustrate the dependence on the initial configurations of the Lloyd iterations and we also address the boundary effects. In Section 5, the qualities of the converged three-dimensional CVT and CVDT are discussed and we give hints on the possible application to high quality meshing. Finally, some conclusions are drawn in Section 6. 2. B A S I C S O F T H E C V T A N D T H E L L O Y D I T E R A T I O N Given a set of input points {zi}~ belonging to a domain l~/ C R k, the Voronoi region corresponding to the point zi consists of all points in 12 that are closer to zi than to any other point in the set. The set {~}~ forms a partition of ~ and is known as a Voronoi tessellation or Voronoi diagram of ft. The points {zi}~' are called generating points or generators. The dual Delaunay triangulation is formed by connecting pairs of generating points which correspond to adjacent Voronoi regions [14,15]. Opt imal CentroidM Voronoi Tessellations 1357 Given a density function p defined on any region Vi, the mass centroids z~ of Vi is defined by . fv~ YP(Y) dy zi fv~ P(Y)dy " A Voronoi tessellation is a centroidal Vovonoi tessellation (CVT) if zi = z*, for all i, i.e., the generators of the Voronoi regions are themselves the mass centroids of those regions. The dual Delaunay triangulation of CVT is referred to as the centroidal Voronoi-Delaunay triangulation (CVDT). When geometric constraints are present, such as the constraints on the boundary vertices, the corresponding definitions of the CVT and CVDT may be generalized to the constrained CVT and constrained CVDT [3,16] and anisotropic CVTs [17]. 2.1. L loyd I t e r a t i o n for C o n s t r a i n e d Cen t ro ida l Voronoi Tesse l la t ion There are several algorithms known for constructing centroidM Voronoi tessellation of a given set [1,18], and they can be categorized into two approaches: the probabilistic and the deterministic approaches. A representative of the probabilistic algorithms is the method by MacQueen using an elegant random sequential sampling and averaging technique [19]. Faster acceleration of the MacQueen's method and its parallel implementation are given in [18]. The deterministic Lloyd algorithm is the obvious iteration between computing Voronoi diagrams and mass centroids [12]. In quantization and clustering literature, one can also find the related h-means and k-means algorithms for the construction of discrete CVTs [11,20]. Here, we choose to work with the Lloyd's method for the construction of three-dimensional CVT and its dual triangulation CVDT mostly because of interests in the energy minimum and the fact that the Lloyd iteration enjoys the energy decreasing property [1]. Given a three-dimensional bounded domain and a sizing function, to construct a constrained three-dimensional CVT (with the generators on the boundary being fixed), the original Lloyd iteration is efficiently implemented through the following procedure [3]. 0 Construct an initial conforming or constrained Delaunay tetrahedronization of a given 3D domain (upon the input of a surface triangular mesh) under a given sizing field H(p) for any given point p in the domain [3]. Here, H(p) can be prespecified or derived from the sizing of the boundary points. Let all the data of the conforming/constrained boundary tetrahedronization be stored. 1 Construct the Voronoi regions for all interior points that are allowed to change their positions, and construct the mass centers of the Voronoi regions with a properly defined density function p derived from the sizing field H(p). Here, the Voronoi regions can be computed directly from the Delaunay tetrahedronization via the duality relationship. 2 Insert the computed mass centers into the conforming boundary tetrahedronization via the constrained Delaunay insertion procedure [3]. 3 Compute the difference D = ~n=l ]IP~ P~mcll 2, where {Pi} is the set of interior points allowed to change, and {Pimc} is the the set of corresponding computed mass center. 4 If D is less than a given tolerance, terminate; otherwise, return to Step 1. In [3], the Lloyd iteration was used for the tetrahedral mesh generation and optimization. Here, it is applied to construct numerous three-dimensional CVTs and their dual CVDTs in order to substantiate the three-dimensional Gersho's conjecture. 2.2. E n e r g y (or Cos t ) P e r Un i t Vo lume and O p t i m a l C V T s V~ ~ Given a density function p, a tessellation V -( i}1 of the domain ~ and a set of points Z -(zi}~ in ~t, we can define the following cost functional: n J:(V,Z) = ~-~ F(V~, z~) , where F(Vi,zi) = [ . p(x)IIx zil[ 2 dx. (1) i=1 1358 Q. DU AND D. WANG This cost functional is also called by other names such as the energy, the error, and the distortion value depending on the type of applications in mind [1]. The standard CVTs along with their generators are critical points of this cost functional (1), as shown in [1]. The Lloyd iteration described earlier enjoys the property that the related functional or cost is monotonically decreasing throughout the iteration [1] thus providing a possible clue for its convergence. More theoretical analysis on the convergence properties of the LIoyd iteration has recently been made in [13]. The energy (or cost) per unit volume (Ep) for a partition or tessellation {V, Z} is then defined by: n 2/~ 7(Y, Z) D(V,Z)k lai:+~/k' n Here, k is the dimension of the space (k = 3 in our paper), [f~l the volume of f / = Ui=l V/. For a given bounded domain f~ together with a specified density function and a fixed number of generator, an optimal CVT is defined as a global minimum of 5r(V, Z), while the optimal centroidal Voronoi tessellation in a given Euclidean-dimensional space (e.g., the two-dimensional space), asymptotically speaking, is defined as the CVT which has the lowest energy per unit volume among all CVTs that cover the whole space (as the number of generators going to infinity). The optimal CVT concept is closely related to the Gersho's conjecture [8], which states that: asymptotically speaking, all cells of the optimal CVT, while forming a tessellation, are congruent to a basic cell which depends on the dimension. This claim is trivially true in one dimension. It has been proved for the two-dimensional case [9] with the basic cell being the two-dimensional regular hexagon. Cersho's conjecture remains open for three and higher dimensions [11]. In [10], it was shown that the body-centered-cubic (BCC) lattice based CVT enjoys the lowest energy per unit volume among all possible lattice based CVTs. The BCC based CVTs has the energy per unit volume valued at 0.07854, with the basic cell given by the truncated octahedron. Other lattices, such as FCC, A15, and Z have the energy values 0.07874, 0.08098 and 0.08666 respectively. Some basic Voronoi ceils for the lattice based CVTs are given in Figure 1. Figure 1. Basic Voronoi cells of the BCC, FCC, A15, and Z configurations. In [21], variants of A15 were also studied, but they were shown to be inferior to the BCC based CVTs. For nonlattice based or general CVTs, it remains unresolved whether the BCC enjoys the lowest energy per unit volume [11]. One question pertains to the possibility of having the optimal CVT made up by a combination of several types of basic cells. In later sections, we design a series of numerical examples for both lattice and nonlattice based CVTs. The computed energy per unit volume and other related properties and statistics substantiate the claim of the three-dimensional Gersho's conjecture: the BCC based CVT enjoys the lowest energy among all three-dimensional CVTs including both lattice and nonlattice CVTs. Thus, asymptotically speaking, the congruent cell of the optimal CVT is the Voronoi cell of the BCC based tessellation, that is, the truncated octahedron. 3. N U M E R I C A L S U B S T A N T I A T I O N OF T H E 3D G E R S H O C O N J E C T U R E To provide convincing numerical evidence for the Gersho conjecture beyond the lattice based argument, in this section, three-dimensional constrained centroidal Voronoi tessellations are conOptimal Centroidal Voronoi Tessellations 1359 structed via the Lloyd iteration for general nonlattice settings. The experiments are designed based on the rationale that, on one hand, the Lloyd iteration always converges to a CVT in practice for any geometry [12,13], and on the other hand, a clear understanding is needed concerning the effect of the geometry boundary, the distribution of initial points and the mesh topological configuration on the final outcome of the Lloyd iteration. Accordingly, this calls for the use of different-initial arrangements in planning our numerical experiments so that we may examine the structure of the converged CVT, and determine whether the result tends to the expected pattern. For a given domain, the initial input of the Lloyd iteration includes a surface triangulation and a distribution of initial field or interior points. The vertices of the surface triangulation and the interior points together make up all the generators of the Voronoi tessellation. Depending on the input geometry and the initial conditions, the optimal CVTs to be computed can be classified into roughly two categories: lattice based CVTs and nonlattice based or general unstructured CVTs. In the former category, all the initial surface triangulations and the interior field points are generated in ways closely related to some kinds of three dimension lattice; while in the latter case, the initial surface triangulation and the interior point distribution are more arbitrary. For all the numerical examples, we compute several useful and indicative statistics of the resulting CVTs. They include: the energy per unit volume (Ep) discussed earlier, the type-ratio, and the quality of the dual Delaunay triangulation (the CVDT). The type-ratios Tr (ratio of the different types of Voronoi cells), in particular, the ratios of polyhedra with 12, 13, 14, 15, and 16 faces, are given as these Voronoi polyhedra typically account for almost all the Voronoi cells in a given CVT and these five types of cells also occupy most of the domain during the Lloyd iteration. The details of the type-ratio reveal much of the geometric structure of the CVT, for example, the type-ratio of a BCC-based CVTs is 0 : 0 : 1 : 0 : 0, which means that all the Voronoi cells are 14-faced polyhedra, (and careful examinations indicate that they are the truncated octahedra). The numerical values given in the paper for the type ratios are not normalized, so that they only provide a relative percentage, for instance, 0.1 : 0.1 : 1.0 : 0.1 : 0.1 would mean the same type ratio as 0.2 : 0.2 : 2.0 : 0.2 : 0.2. The element quality data on the centroidal Voronoi Delaunay triangulation (CVDT) show that the optimal CVTs in three-dimensional space can lead to high quality tetrahedral meshes. Here, for each tetrahedral element e, we use the quality measure defined by Q(e) = 2x/6R~n/hm~x with Rin being the radius of the inscribed-sphere of the tetrahedron, and hm~× the longest edge length. We start our investigation with simple lattice based CVTs and then move onto more general CVTs. We note that the statistics are collected via the implementation of algorithms specifically designed for our numerical investigations. 3.1. P r o p e r t i e s of La t t i c e s B a s e d C V T s Lattices based CVTs were widely used in vector quantization, coding theory, bubble and foam geometry [10,21-24]. Various methods for the construction of different kinds of lattices and their quantizations or CVTs have been reported in [21,23,25,26]. In [10], it was shown that the optimal lattice quantizer in three dimensions is the body-centered cubic (BCC) lattice based quantizer, i.e., the optimal lattice-based CVT cell is given by the truncated octahedron. Recently, a variant of the A15 lattice (so-called Weaire-Phelan partition) [23,25,27] was proposed in connection with the Kelvin conjecture (on the partition of R 3 with the least surface area). It has been shown in [10,23,25] that this A15 variant is inferior to the BCC in terms of both the value of Ep and the normalized moment of inertia (NMI). To provide further numerical substantiation to the Gersho's conjecture, we design the following experiments: first, we compute the optimal CVTs from a set of initial conditions for a cubic domain via the Lloyd iteration. The cubic domain is scaled to be [0, 10] ~ without loss of generality. We first construct a series of lattice-based CVTs including cells represented by the FCC (the face centered Cubic lattice), BCC, A15, and Z configurations. The A15 configuration here refers to the original partition, not the Weaire-Phelan partition. As most of the lattice-based CVTs are 1360 Q. Du AND D. WANG (a) (b) (c) (d) Figure 2. FCC lattice: (a) the surface triangulation of the original cube [0, 1013; (b) the FCC based surface triangulation for elements inside a smaller cube [2, 8]3; (c) a cutt ing view of the interior Delaunay triangulation of the FCC lattice points; (d) a cutt ing view of the FCC based CVDT. not exactly cube-filling, to minimize the boundary effect, we only collect the needed statistics from a smaller interior region away from the boundary. We also carry out similar computations for domains other than a cube such as cylinders and other composed domains. With the cubic domain as a primary example, our computational procedure involves the following main steps: (1) the surface triangulation of the domain is first constructed with uniform sizing (as shown in Figure 2a for the cubic domain [0, 1013); (2) then, for a lattice based configuration, the inner field points are generated using the corresponding lattice (with the same sizing as the surface triangulation); (3) with the inner field points obtained above, an initial constrained three-dimensional Delaunay triangulation is generated via the method of constrained Delaunay insertion and associated boundary recovery procedure [3,28,29]; (4) the elements contained in a smaller interior domain (for example, the smaller cube [2, 8] 3 inside the original cubic domain [0, 10] 3) are gathered to form a Delaunay triangulation whose dual configuration gives a Voronoi tessellation corresponding to the pure lattice points, without using any boundary points of the original larger domain (take the FCC as an example, Figure 2b shows a surface triangulation of the FCC based Delaunay mesh, while one of its cutting view is shown in Figure 2c); (5) finally, starting from this Delaunay mesh, the Lloyd iteration is performed until convergence to obtain the final CVT and its dual CVDT ~Figure 2d shows a cutting view of an example given by the FCC-based CVDT). We note that in our implementation of the Lloyd iteration, no movement of the boundary generators is allowed in order to preserve the boundary integrity. The use of the lattice based initial field distribution is to study the stability of corresponding lattices, i.e., whether a lattice-based Voronoi tessellation leads asymptotically to a local or global minimum of the energy (and thus Ep), as well as their sensitivity to the boundary surface triangulation. FCC. The face-centered-cubic (FCC) lattice is formed by the vertices of the regular cubic ceils and their face centers with the corresponding Voronoi cell being a rhombic dodecahedron and the type-ratio 1 : 0 : 0 : 0 : 0. The Voronoi tessellation of a FCC lattice is a CVT with Ep -0.07874 and its dual Delaunay triangulation has two kinds of basic elements with element qualities at Q = 0.656 and Q = 1.00 respectively. The dihedral angles corresponding to the two basic elements are [4(54.735), 1(90), 1(109.47)] and [6(70.528)]. The average mesh quality number is 0.7823. We note that in [30], it has been argued that the FCC was a preferred choice, comparing to the BCC, for an optimal finite-element mesh generation in terms of better approximation error bounds. To study the stability of the FCC lattice based CVT, we first either make some random perturbations to the interior lattice points that lead to an initial FCC lattice Voronoi tessellation, or Optimal Centroidal Voronoi Tessellations 1361 replace these points with some other point distributions, then apply the Lloyd iteration to construct constrained CVTs. Similar experiments have been reported in [1]in the two-dimensionai cases where it was demonstrated that the regular square lattice is less stable and it often transforms to the hexagonM lattice under initial perturbations and the Lloyd iteration. For our experiments, we denote the case where the initial perturbations move the vertices by less than 10% of the local sizing as PS and the case where the perturbations move the vertices to as much as 50% of the local sizing as PL. The case PC denotes the case where the interior points are replaced by small random perturbations to points on a Cartesian grid with comparable number of vertices through an appropriate scaling. Case PS PL PC Table 1. tions. Statistics of FCC based CVTs and CVDTs with different initial distribu-